subroutine output
use cell
use iounits
use general
use wallinfo
use waterinfo
use H2MolInfo
implicit none
   !----------------------------------------------------------------------------
   integer :: i, II, im, iatom, idummy, imw, imwmax, ibdw = 0, wid = 0
   real(q) :: xlo, xhi, ylo, yhi, zlo, zhi
   !----------------------------------------------------------------------------
   xlo = 0.D0
   xhi = box(1)
   ylo = 0.D0
   yhi = box(2)
   zlo = 0.D0
   zhi = box(3)
   !
   open ( ioout, file = outfile, action = 'write', status = 'unknown' )
   write( ioout, 10 ) trim(title)
   ! headers
   ibdw = 0
   idummy = natom
   if ( LAddWalls ) idummy = idummy + sum(WallRemain)
   if ( lwater    ) idummy = idummy + nH2Oatom
   if ( LH2Mol    ) idummy = idummy + nH2atom
   write( ioout, 20 ) idummy, "atoms"
   idummy = nbonds
   if ( lwater    ) idummy = idummy + nH2Obond
   if ( LH2Mol    ) idummy = idummy + nH2Bond
   write( ioout, 20 ) idummy, "bonds"
   idummy = nangls
   if ( lwater    ) idummy = idummy + nH2Oangle
   write( ioout, 20 ) idummy, "angles"
   write( ioout, 20 ) ndihrs, "dihedrals"
   write( ioout, *  )
   idummy = natmtyp
   if ( LAddWalls ) idummy = idummy + size(WallIndex, dim=1)
   if ( lwater    ) idummy = idummy + 2
   if ( LH2Mol    ) idummy = idummy + 1
   write( ioout, 20 ) idummy,  "atom types"
   idummy = nbndtyp
   if ( lwater    ) idummy = idummy + 1
   if ( LH2Mol    ) idummy = idummy + 1
   write( ioout, 20 ) idummy, "bond types"
   idummy = nangtyp
   if ( lwater    ) idummy = idummy + 1
   write( ioout, 20 ) idummy, "angle types"
   write( ioout, 20 ) ndihtyp, "dihedral types"
   write( ioout, *  ) 
   write( ioout, 30 ) xlo, xhi, "xlo", "xhi"
   write( ioout, 30 ) ylo, yhi, "ylo", "yhi"
   write( ioout, 30 ) zlo, zhi, "zlo", "zhi"
   write( ioout, *  )
   !
   ! atom property section
   ! mass
   write( ioout, 40 ) "Masses"
   do i = 1, natmtyp
      write( ioout, 50 ) i, EMass(i)
   enddo
   idummy = natmtyp
   if ( LAddWalls ) then
      do i = 1, size(WallIndex, dim=1)
         idummy = idummy + 1
         write( ioout, 50 ) idummy, WallMass(WallIndex(i))
      enddo
   endif
   if ( lwater    ) then
      write( ioout, 50) idummy+1, H2OMass(1)
      write( ioout, 50) idummy+2, H2OMass(2)
      idummy = idummy + 2
   endif
   if ( LH2Mol    ) then
      write( ioout, 50) idummy+1, H2Mass(1)
   endif
   !
   ! position: atom
   write( ioout, 40 ) "Atoms"
   do i = 1, natom
      write( ioout, 60 ) i, molID(i), attyp(i), atchg(i), atpos(:,i)
   enddo
   iatom = natom
   im    = maxval(molID)
   !
   idummy = natmtyp
   wid    = 0
   if ( LAddWalls ) then
      imwmax = 0
      do i = 1, nWallAtom
         if ( WallRemain(i).eq.1 ) then
            iatom = iatom + 1
            write( ioout, 60 ) iatom, im+WallMolID(i), idummy+WallTyp(i), WallCHG(WallIndex(WallTyp(i))), WallPos(:,i)
            imwmax = max(imw, imwmax)
         endif
      enddo
      idummy = idummy + maxval(WallTyp)
      im = im + maxval(WallMolID)
      wid = im
   endif
   if ( lwater ) then
      do i = 1, nH2Oatom
         iatom = iatom + 1
         write(ioout, 60) iatom, im+waterMolID(i), idummy+watertype(i), waterchg(i), waterpos(:,i) 
      enddo
      idummy = idummy + maxval(watertype)
      im = im + maxval(waterMolID)
   endif
   if ( LH2Mol ) then
      do i = 1, nH2atom
         iatom = iatom + 1
         write(ioout, 60) iatom, im+H2MolID(i), idummy+H2MolType(i), H2MolChg(i), H2MolPos(:,i) 
      enddo
      idummy = idummy + maxval(H2MolType)
      im = im + maxval(H2MolID)
   endif
   ! molecular topology section
   ! bonds
   write( ioout, 40 ) "Bonds"
   do i = 1, nbonds
      write( ioout, 70 ) i, bonds(:,i)
   enddo
   idummy = natom
   iatom  = natom
   if (LAddWalls) idummy = idummy + sum(WallRemain)
   if ( lwater ) then
      do i = 1, nH2Obond
         write(ioout, 70) nbonds+i, nbndtyp+1, waterbond(:,i)+idummy
      enddo
      ibdw = 1
      nbonds = nbonds + nH2Obond
      iatom  = iatom  + nH2Oatom
   endif
   if ( LH2Mol ) then
      do i = 1, nH2Bond
         write(ioout, 70) nbonds+i, nbndtyp+ibdw+1, H2MolBond(:,i) + iatom
      enddo
   endif
   ! angles
   write( ioout, 40 ) "Angles"
   do i = 1, nangls
      write( ioout, 80 ) i, angls(:,i)
   enddo
   if ( lwater ) then
      do i = 1, nH2Oangle
         write(ioout, 80) nangls+i, nangtyp+1, waterangle(:,i) + idummy
      enddo
   endif
   ! dihedrals
   write( ioout, 40 ) "Dihedrals"
   do i = 1, ndihrs
      write( ioout, 90 ) i, dihrs(:,i)
   enddo
   !
   ! Force field section
   ! bond coefficients
   write( ioout, 40 ) "Bond Coeffs"
   do i = 1, nbndtyp
      write( ioout, 95 ) i, bondcoef(:,i)
   enddo
   if (lwater) write( ioout, 95 ) nbndtyp+1, HOBondCoeff
   if (LH2Mol) write( ioout, 95 ) nbndtyp+ibdw+1, HHBondCoeff
   ! angle coefficients
   write( ioout, 40 ) "Angle Coeffs"
   do i = 1, nangtyp
      write( ioout, 95 ) i, anglcoef(:,i)
   enddo
   if (lwater) write( ioout, 95 ) nangtyp+1, HOHAngCoeff
   ! dihedral coefficients
   write( ioout, 40 ) "Dihedral Coeffs"
   do i = 1, ndihtyp
      write( ioout, 99 ) i, torscoef(:,i)
   enddo
   ! pair coefficients
   write( ioout, 40 ) "Pair Coeffs"
   do i = 1, natmtyp
      write( ioout, 95 ) i, epison(i), sigma(i)
   enddo
   idummy = natmtyp
   if ( LAddWalls ) then
      do i = 1, size(WallIndex, dim=1)
         idummy = idummy + 1
         write( ioout, 95 ) idummy, WallLJPM(:,WallIndex(i) )
      enddo
   endif
   ibdw = 0
   if ( lwater ) then
      write( ioout, 95 ) idummy+1, H2OPair(:,1)
      write( ioout, 95 ) idummy+2, H2OPair(:,2)
      ibdw = 2
   endif
   if ( LH2Mol ) then
      write( ioout, 95 ) idummy+ibdw+1, HHPair(:,1)
   endif
   close( ioout )
   if ( LAddWalls ) then
      write( stdout, 100 )
      write( stdout, 110 )
      write( stdout, 120 ) "boundary p p f #p p fs"
      write( stdout, 130 ) "group upper molecule",  wid-1
      write( stdout, 130 ) "group lower molecule",  wid
      write( stdout, 140 ) "group liquid molecule", wid-2
      write( stdout, 150 ) "group walls union upper lower"
      write( stdout, 150 ) "fix freezewalls walls setforce 0.0 0.0 0.0"
      write( stdout, 125 ) "neigh_modify check yes exclude type ",  natmtyp+1, natmtyp+1
      write( stdout, 150 ) "kspace_modify   slab 3.0"
      write( stdout, 120 ) "compute  MyTemp liquid temp"
      write( stdout, 150 ) "compute  MyPres all pressure MyTemp"
      write( stdout, 150 ) "thermo_modify  temp MyTemp press MyPres"
      write( stdout, 150 ) "#pair_modify  tail yes SHOULD NOT be used!"
      write( stdout, 100 )
   endif
   if (lwater) then
      write( stdout, 135 ) "group water type", idummy+1, idummy+2
   endif
   if (LH2Mol) then
      write( stdout, 135 ) "group h2mol type", idummy+ibdw+1
   endif
   !
10 format(A, /)
20 format(I20,2x,A)
30 format(F16.8,1X,F16.8,2x,A3,1x,A3)
40 format(/,A,/)
50 format(2x,I3,2x,F12.6)
60 format(1x,I7,1x,I6,1x,I4,1x,F15.8,3(1x,F15.8))
70 format(1x,I7,1x,I5,2(1x,I7))
80 format(1x,I7,1x,I5,3(1x,I7))
90 format(1x,I7,1x,I5,4(1x,I7))
95 format(1x,I5,2(1x,F12.6))
99 format(1x,I5,5(1x,F12.6))
100 format(/, 10x,80("="))
110 format( 10x, "Additional information needed in the run script:", / )
120 format( 10x, A, / )
125 format( 10x, A, 2(1x,I3),/)
130 format( 10x, A, 2x, I5 )
135 format( 10x, A, 2(2x,I5))
140 format( 10x, A, 2x, "<=", I4 )
150 format( 10x, A     )
return
end subroutine
